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ABSTRACT 


Previous  experimental  investigations  have  shown  that  the  characteristics  of  flow 
about  a  circular  cylinder  immersed  in  a  time-dependent  flow  exhibit  cycle-to-cycle 
variations.  These  variations  have  been  attributed  to  the  variations  in  the  spanwise 
coherence,  aspect  ratio,  nonuniformity  of  the  flow,  and  random  disturbances  in  the 
ambient  flow.  A  theoretical  investigation  was  undertaken  to  examine  the  stability  of 
the  flow  characteristics  in  terms  of  the  initial  state  of  the  vortices.  An  idealized  model 
has  been  devised  and  the  position  of  the  vortex  was  varied  systematically.  The  results 
have  shown  that  finite-precision  information  about  the  characteristics  of  the  flow  does 
not  lead  to  finite-precision  information  at  a  later  stage.  In  fact  the  advection  of  the 
vortices  can  give  rise  to  chaotic  behavior  in  the  calculated  lift  and  drag  forces  and  in 
the  velocity  field.  It  is  concluded  that  the  cycle-to-cycle  variations  are  not  entirely  due 
to  lack  of  spanwise  coherence  and  that  they  are  mostly  a  consequence  of  the  chaotic 
motion  which  can  result  from  the  advection  of  the  vortices  in  a  time-dependent  flow. 
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I.  INTRODUCTION 


In  the  past  ten  years  there  has  been  much  effort  to  analyze  the  separated  time 
dependent  How  about  bluff  bodies  in  general  and  circular  cylinders  in  particular. 
Shortfalls  in  the  M orison  Equation  to  adequately  predict  in-line  forces  as  weil  as  the 
advent  of  sophisticated  computational  and  experimental  techniques  have  stimulated  a 
great  deal  of  research  in  this  area.  The  accumulation  of  this  numerical,  experimental 
and  visual  data  has  generated  the  need  for  a  deeper  understanding  of  the  motion  of 
vortices  in  time  dependent  How.  Specifically  in  oscillating  flow  the  stability  of  the 
vortices  generated  is  of  particular  importance.  The  significance  of  the  stability  of  the 
vortices  is  a  result  of  its  possible  exposition  of  documented  cycle-to-cycle  variations  in 
measured  forces  on  the  cylinder. 

In  this  application,  stability  is  not  a  characteristic  of  the  vortex  life  or  demise  but 
of  the  vortex  motion  or  path.  A  system  is  referred  to  as  stable  in  any  particular  state 
of  equilibrium  if.  after  being  given  a  finite  perturbation,  it  tends  to  return  to  the  state 
of  equilibrium  existing  before  it  was  disturbed.  On  the  contrary  an  unstable  system 
subjected  to  a  finite  perturbation  catastrophically  fails  to  return  to  its  initial  state  of 
equilibrium.  Many  real  systems  in  nature  do  not  fall  into  either  of  these  two  categories 
though.  In  these  systems  there  are  certain  times  when  finite  disturbances  do  not  lead 
to  finite  results  nor  catastrophic  or  runaway  reactions.  These  complicated  systems  are 
described  as  neither  stable  nor  unstable  but  as  systems  susceptible  to  chaos. 

Whether  or  not  the  stability  of  the  vortex  motion  in  oscillating  flows  can  be 
categorized  as  chaotic  requires  not  only  an  analysis  of  the  vortex  motion  but  also  a 
dear  definition  and  understanding  of  chaos.  Unfortunately  the  pursuit  of  this 
understanding  remains  the  subject  of  ongoing  research  and  a  concise  definition  is  yet 
elusive.  Strides  have  been  made  however  in  identifying  chaotic  behavior.  Systems 
studied  that  demonstrate  chaotic  behavior  could  best  be  described  as  bifurcating 
systems  that  operate  along  a  razor's  edge  or  precipice.  The  slightest  disturbance  will 
pusii  the  system  one  way  or  the  other  and  will  lead  to  decidedly  different  (but  not 
unstable)  results.  Additionally,  the  systems  studied  reveal  that  the  state  of  chaotic 
behavior  is  no:  all  encompassing  but  exists  more  in  regions  of  susceptibility  to  chaotic 
behavior  preceded  and  followed  by  regions  of  stable  behavior. 


1  ! 


Th.e  application  of  this  analysis  of  chaotic  behavior  to  vortex  motion  in 
oscillating  flow  is  particularly  useful  due  to  the  complexity  of  the  factors  driving  the 
vortex  trajectory.  The  coupling  of  relative  position  to  other  vortices  (real  or  image)  in 
-a  time  dependent  How  is  both  intricate  and  inseparable.  Additionally,  the  flow  field 
about  a  bluff  body  near  the  surface,  combined  with  the  problems  of  multiple  vortices, 
viscous  dissipation,  wall  reflection  and  boundary  layer  and  limited  coherence  length,  is 
extremely  intricate.  Consequently,  the  motion  of  the  vortex  at  times  becomes  subject 
to  finite  perturbations  without  finite  results.  The  specific  coupling  that  makes  the 
system  susceptible  to  chaotic  behavior  gives  rise  to  the  regions  of  chaos  preceded  and 
followed  by  regions  of  stability. 

In  order  to  understand  these  complicated  causes  and  effects  first  an 
understanding  of  the  motion  of  a  single  potential  vortex  in  oscillating  flow  about  a 
bluff  body  is  required.  Therein  lies  the  impetus  of  this  study.  In  an  ideal  two 
dimensional  potential  flow,  wall  refection,  wail  boundary  layers  and  coherence  length 
are  eliminated.  Additional  vortices  can  be  added  to  study  mutual  induction  effects. 
Further,  advanced  models  of  the  vortex  (i.e.  Rankine  and  Gaussian  time  dependent 
models)  can  clarify  the  broad  effects  of  viscous  dissipation. 

The  scope  of  this  study  includes  the  potential  flow  of  one  and  two  vortices  placed 
in  the  vicinity  of  a  circular  cylinder  in  oscillating  flow.  Emphasis  is  placed  on 
identifying  the  onset  of  chaotic  behavior  and  its  effects  on  the  flow  characteristics. 


II.  BACKGROUND 


The  concept  of  chaos  in  dynamic  systems  is  by  no  means  new.  However,  the 
detailed  investigation  as  to  the  onset  and  mechanisms  of  chaos  is  still  comparatively 
young.  The  primary  difficulty  in  the  investigation  of  chaos  is  the  fact  that  by  definition 
chaotic  behavior  is  non-integrable.  Because  of  this  the  first  attempts  to  define  the 
concept  of  chaos  were  undertaken  by  mathematicians.  In  1984  Holmes  [Ref.  1] 
investigated  chaotic  motion  in  forced  oscillations.  Here  the  onset  of  chaos  was  linked 
to  periodicity  or  aperiodic  behavior  of  the  function  through  scrutiny  of  the  harmonics. 
Although  this  is  a  credible  definition  its  application  appears  limited  to  the  idealized 
mathematical  system  from  which  it  was  developed.  In  1980  Aref  [Ref.  2]  expanded  on 
this  definition  and  the  use  of  Poincare'  maps  by  illustrating  chaotic  behavior  in  paths 
of  finite,  ideal,  and  controllable  vortices.  His  motivation  for  the  investigation  was  the 
onset  and  mechanisms  of  chaos.  Therefore,  while  the  idealized  system  that  he  designed 
was  useful  for  demonstrating  chaotic  behavior,  its  application  to  real  world  dynamics 
was  severly  limited. 

Other  efforts  were  made  to  try  to  link  the  study  of  chaos  to  systems  more  closely 
resembling  those  found  in  nature.  In  1984  Hardin  &  Mason  [Ref.  3]  investigated  the 
two  dimensional  system  of  axial  vortices  in  a  pipe.  This  was  the  first  attempt  to 
examine  a  system  that  resembled  real  phenomena  and  that  also  underwent  chaotic 
behavior.  All  previous  studies  examined  chaotic  behavior  by  generating  systems 
susceptible  to  chaos  and  then  extrapolating  the  results  to  possible  real  world  problems. 

Finally,  in  1986  Dowell  &  Pezeshki  [Ref.  4|  investigated  the  one  dimensional 
Duffings  Equation  and  the  onset  of  chaos.  In  its  relation  to  real  systems  Duffings 
Equations  could  be  used  to  describe  the  snap  buckling  of  a  beam  and  other 
phenomena.  In  the  beam  the  onset  of  chaos  is  easily  observed  but  the  integration  of 
the  equation  to  predict  or  study  the  mechanisms  of  chaos  is  somewhat  abstract.  Of 
particular  interest  in  this  study  was  the  examination  of  the  onset  of  chaos  by  Phase 
Plane  Plots  (displacement  -  velocity  plots). 

By  reasons  of  their  treatment  of  fluid  dynamics  in  general  and  vortex  motion  in 
particular.  Aref s  chaotic  advection  and  Hardin  &  Mason's  vortices  in  a  pipe  have 
direct  application  to  this  investigation  and  will  be  discussed  in  greater  detail. 


A.  CHAOTIC  ADVECTION 

Aref  [Ref.  2]  devised  a  numerical  model  to  represent  the  motion  of  an 
incompressible  inviscid  fluid  within  a  circular  bounding  contour.  A  stirring  agitator 
was  modelled  as  a  potential  point  vortex.  The  agitator(s)  provided  the  potential  flow 
within  the  boundary.  A  global  Eulerian  view  of  an  array  of  markers  was  traced  to 
examine  and  verify  the  poor  efficiency  of  the  laminar  mixing. 

In  a  follow-up  study,  he  examined  the  case  of  two  fixed  agitators  (point  vortices) 
in  the  fluid  that  generated  "piecewise-constant  stirrer  motion."  The  vorticity  of  the 
agitators  was  controllable  (on  or  off)  and  switched  back  and  forth  at  prescribed  times 
to  generate  unsteady  flow.  Again  a  marker  array  of  particles  was  traced.  The  mixing 
was  clearly  superior  and  Aref  classified  the  motion  as  chaotic.  The  parameters  that 
dictate  the  switching  of  the  vortices  and  therefore  the  chaotic  regions  were  cited  as 
mechanisms  for  chaos. 

Extrapolations  to  real  phenomena  include  enhanced  mixing  chambers  with 
controllable  agitators  and  similarities  to  gravitational  systems. 

B.  VORTICES  IN  A  PIPE 

The  motion  of  axial  vortices  was  investigated  by  Hardin  and  Mason  [Ref.  3] 
because  of  its  widespread  appearance  in  physical  situations.  These  vortices  occur  in 
virtually  every  flow  that  changes  axial  directions  (i.e.  rivers,  pipes,  lung  bronchial  tubes, 
etc.).  In  these  physical  systems  pairs  of  counterrotating  vortices  are  produced.  This 
study  investigated  the  effects  of  one  and  two  pairs  of  equal  strength  counterrotating 
two  dimensional  potential  vortices  enclosed  in  a  circular  boundary.  Vortex  advection 
by  induced  velocities  (real  or  image)  was  determined  and  displayed  by  means  of 
Eulerian  paths  of  the  vortices. 

Controlling  factors  on  this  investigation  were  initial  locations  of  the  vortices 
within  the  boundary.  For  the  case  of  single  pair  of  vortices  the  paths  were  found  to  be 
periodic  or  stable.  In  the  case  of  two-vortex  pairs  (four  vortices  within  the  boundary), 
the  periodicity  of  the  vortex  path  and  flow  was  a  function  of  where  the  vortices  were 
placed  initially.  This  supports  the  concept  of  chaotic  motions  occurring  within  regions 
of  susceptibility.  It  was  found  that  within  these  regions  of  susceptibility  that  small 
perturbations  (0.001  radii  of  the  enclosing  boundary)  were  sufficient  to  drive  the  vortex 
path  aperiodic  and  demonstrate  chaotic  behavior.  Hardin  &  Mason  [Ref.  3]  do  note 
that  experimental  evidence  indicates  that  bifurcating  or  branching  flows  are 
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quantitatively  more  stable  in  real  viscous  flows  than  these  numerical  results  suggest. 
This  is  attributed  to  the  complexity  of  this  secondary  flow  as  well  as  viscous 
dissipation,  vortex  decay  and  coring  influences.  It  does  not,  however,  diminish  the 
validity  of  the  study  or  the  qualitative  results  of  the  susceptibility  of  fluid  flows  to  the 
onset  of  chaos. 
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III.  ANALYSIS 

A.  INTRODUCTION 

In  1976  Sarpkaya  [Ref.  5]  presented  extensive  data  on  oscillating  flow  about 
circular  cylinders.  In  addition  to  modifying  the  long  standing  Morison  Equation,  this 
data  also  provided  extreme  insight  into  transverse  forces  on  the  cylinder.  Performed  at 
Reynolds  numbers  of  prototype  magnitude  and  across  a  wide  spectrum  of  Keulegan- 
Carpenter  [Ref.  6)  numbers,  these  experiments  provided  a  data  base  to  which  all 
subsequent  publications  and  literature  on  the  subject  has  been  based.  Of  particular 
significance  and  not  predicted  by  the  Morison  Equation  was  the  documentation  of 
cycle-to-cycle  variations  in  pressures  and  forces  across  the  cylinder  in  pure  harmonic 
flow  (see  Figure  E.l).  These  variations,  particularly  significant  in  the  transverse 
direction,  were  attributed  primarily  to  limited  coherence  lengths  and  to  a  lesser  extent 
wall  reflections  and  boundary  layers. 

The  distorting  effects  of  coherence  length  on  the  wake  of  a  circular  cylinder  in 
uniform  flow  was  clearly  discussed  by  Gerlach  [Ref.  7]  and  Graham  [Ref.  8.J  In 
oscillating  flows,  however,  the  von  Karman  street  forms  perpendicular  or  transverse  to 
the  flow  direction  in  the  range  8  <  Kc  <13.  This  makes  experimentation  to  identify  the 
coherence  length  extremely  difficult.  In  an  effort  to  specifically  examine  the  effects  of 
spanwise  coherence,  O'Keefe  [Ref.  9]  in  1986  undertook  similar  experiments  to 
Sarpkaya  [Ref.  5]  but  with  the  minimum  aspect  ratio  (L.'D  =  2.0)  that  could  be 
accurately  instrumented.  These  nearly  two  dimensional  experiments  provided 
extremely  accurate  data  that  was  not  blurred  by  the  influence  of  a  limited  coherence 
length.  Cycle-to-cycle  variations  in  in-line  and  transverse  forces  remained  though; 
particularly  in  the  Keulegan-Carpenter  number  range  8  <  K  <13. 

Knowing  that  at  any  given  point  in  time  the  forces  acting  on  the  cylinder  is 
governed  primarily  by  one  dominant  vortex,  a  detailed  investigation  of  one  potential 
vortex  in  the  vicinity  of  a  circular  cylinder  in  harmonically  oscillating  flow  was  dictated. 
After  an  investigation  of  the  single  vortex  case  the  follow-up  was  to  investigate  two 
potential  vortices  in  the  vicinity  of  a  circular  cylinder. 
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B.  POTENTIAL  FLOW  MODELING 
1.  One  Vortex  Model 

In  the  complex  plane  the  potential  flow  field  for  a  vortex  in  the  vicinity  of  a 
circular  cylinder  in  harmonically  oscillating  flow  is  given  by  the  following  expression: 

ff z)  =  U  vsin(cot)(z  +  — )  +  i— ln(z-z  )-i— ln(z-^-)  (3.1) 

max  c  2n  0  z 

O 

where:  c  -  radius  of  the  cylinder 

T  -  strength  of  the  vortex 
(0  -  frequency  of  the  oscillating  flow 
L'max  -  maximum  amplitude  of  the  flow  velocity 
z  -  location  of  the  vortex 

O 

In  nondimensional  form  this  expression  becomes: 


f(Q=sin( 2k t){ £  +  4-)  +  iKln( £  -  C  )  -  iKln{ q  -  ^-) 
■»  Co 
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(3.2) 


where:-  K  =  r/(2TCcUmax) 

r  =  t/T  (T  -  period  of  the  flow) 

C  =  z/c 

This  expression  can  be  differentiated  with  respect  to  q  to  determine  the  flow 
field  velocities.  At  C  =  C0  the  expression  for  the  velocity  of  the  vortex  becomes: 


u  —  iv  =  sin( 2rcr )( l  - 


iK 


(3.3) 


where  r  is  the  radius  between  the  real  and  image  vortex  and  is  given  bv: 


(Note  that  the  second  term  of  eqn  3.2  is  dropped  since  a  vortex  does  not 
impart  velocity  to  itself.) 

Calculations  for  the  forces  acting  on  the  cylinder  were  made  by  integrating  the 
Blasius  Force  Fquation: 


D  -  iL  =  — 1- ipf(— )2dz  +  ip— ^fdz 
2  7  dz  'dr 


(3.5) 


Nondimensionalized  and  evaluated  at  £j*>*£  ,  this  equation  reduces  to  the 
following  form: 


D'  -  iL'  =  -2^—  cos(2ttt)  +  (£  "r-) 
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(3.6) 
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where:  K  =  (l'  yT)'(2c) 


2.  Two  Vortex  Model 

Invoking  the  same  procedure  as  the  one  vortex  model  it  can  be  shown  that  the 
velocity  of  two  vortices  A  and  B  in  the  vicinity  of  a  the  circular  cylinder  is  given  by: 


&?\  -  u-iv  *  sin(2nt)(l--4)-(- — ^^) 


iK— )  +  (— ■  (3.7) 
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ttI  =  U-iv  =  sin(  2tx  r )( l  —75)  —  (3 - ^ - )  ~  (  -  irx; — )  +  (- - — )  (3.8) 
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3.  Rankine  Vortex  Model 

A  standard  rankine  vortex  model  was  used  with  linear  velocity  distribution 
between  0  and  rc0rc-  l  he  value  oi  tltc  cue  radius  was  droved  to  be  n  lu  times  tire 
radius  of  the  cylinder  generating  the  sortu.es 
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4.  Gaussian  Time  Dependent  Vortex  Model 

A  standard  Gaussian  time  dependent  vortex  model  was  used  where  the 
generalized  expression  for  the  velocity  distribution  is: 


_  _ r 


2ttr 


Here  t  is  not  the  clock  time  of  the  How  but  the  life  time  of  the  vortex.  It  has 
been  shown  that  differentiating  this  with  respect  to  the  radius  yields  r  =  2.24 v  ( vtv >  as 
the  expression  for  the  core  radius  of  the  vortex.  Nondimensionalizing  eqn  3.S  in  terms 
of  the  oscillating  flow  parameters,  and  examing  the  velocity  of  a  vortex  from  eqn  3.3. 
the  velocity  of  a  Gaussian  time  dependent  vortex  is  given  by: 


u~ 


i  ;l-  r  "v  (K  Re  ) 

iv  =  sin(2rcr)<  1  ~-rj~  1 1  -  expi  -  - 

q2  rv  16(tv  T) 


.)}  (3.10) 


For  the  modeling  the  Reynolds  Number  was  chosen  to  be  Re=  10.000  and  tv 
was  set  such  that  the  core  of  the  vortex  was  r  =0.10  at  inception. 

5.  Numerical  Integration  and  Accuracy 

In  this  analysis  a  finite  differencing  technique  was  employed.  The  study  of 
potentially  chaotic  systems  requires  a  high  degree  of  accuracy.  Due  to  its  guaranteed 
stability  and  overall  efficiency,  the  Adams-Moulton  predictor  corrector  method  was 
selected  to  convect  the  vortices: 


znp  =  zn-l  +  <55  Vl  ”  59  +  37  %-3  “  9  %-4) 


7n  =  zn-l  '  -9 


(3.11 ) 


(3.12, 


In  addition  to  ensured  stability,  tins  method  provided  fourth  order  accuracy 
with  local  error  Oh5).  With  1440  time  steps  per  cycle,  the  local  error  remained  at  the 
limits  of  double  precession  machine  accuracy.  Initial  startup  was  made  using  the  lour*;; 
order  Rimge-Kutta  scheme.  All  continuations  were  based  on  the  last  four  calcul.red 
locations  of  the  vortices,  therebv  ensuring  accuracv  and  duplication  at  am  cede. 
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While  discussing  accuracy  it  must  be  noted,  as  mentioned  earlier,  that  by 
definition  chaotic  behavior  is  non-integrable.  Whatever  scheme  used  can  only  illustrate 
qualitatively  the  onset  of  chaos.  If  used  for  comparison  universally,  all  methods  yield 
equally  valid  qualitative  results.  The  effort  devoted  here  to  accuracy  is  therefore  to 
ensure  the  validity  of  the  results  and  add  to  overall  integrity  of  this  report. 


IV.  RESULTS 


A.  GENERAL 

From  Sarpkava  [Ref.  5]  and  other  studies,  it  is  concluded  that  the  region  of 
special  interest  in  the  field  of  oscillating  Hows  is  confined  to  8<  Kc<  13,  an  interval  in 
which  a  half  transverse  vortex  street  is  generated.  It  is  because  of  this  reason  that  in 
the  present  numerical  experimentation  Kc  was  chosen  equal  to  10.0. 

The  vortex  path  is  a  function  of  both  the  imposed  ambient  flow  field  and  the 
velocities  induced  by  the  vortices  (real  or  image).  This  coupling  was  first  investigated 
through  the  use  of  a  single  vortex  in  uniform  How  in  the  vicinity  of  a  circular  cylinder. 
Figure  E.2  illustrates  the  path  of  the  vortex  for  various  vortex  strengths  (K  =  0.5,  0.625, 
0.75).  It  is  seen  that  the  capture  by  the  image  vortex  (partial  looping)  lasts  only  until 
the  How  field  around  the  cylinder  pushes  the  vortex  away  from  the  image.  It  is  also 
seen  that  this  capture  is  a  function  of  the  initial  location  of  the  vortex.  When  the 
vortex  is  captured,  sharp  cusps  appear  where  the  direction  of  motion  of  the  vortex  is 
quickly  reversed.  These  characteristics  played  even  more  decisive  roles  in  the  unsteady- 
case.  With  the  exception  of  the  foregoing,  the  vortex  strengths  used  in  the  numerical 
experimentation  were  all  set  at  K  =  0.50. 

B.  ONE  VORTEX  MODEL 

Initial  experimentation  included  varying  the  starting  locations  and  times  of  the 
single  point  vortex  in  order  to  observe  the  vortex  path  and  force  coefficients.  Figures 
E.3  through  E.6  show  the  force  coefficients  for  various  starting  locations  and  times  and 
varying  numbers  of  cycles.  These  figures  clearly  show  cycle-to-cycle  variations  in  both 
the  inline  and  transverse  force  coefficients.  In  comparison  to  experimental  force  traces 
shown  in  Figure  F.l.  of  comparable  Kc.  the  variations  are  of  similar  magnitude  and 
shape. 

Understanding  the  origin  of  these  cycle-to-cycle  variations  required  a  scrutiny  of 
the  vortex  paths  themselves.  Figures  E.7  and  E.8  show  two  starting  positions  of  the 
vortex  of  (0,  1.5083)  and  (0,  1.5084).  The  oscillating  flow  velocity  was  set  at  U  =0  + 
(superscript  indicates  the- initial  x-dircction).  Figures  E.7(a)  and  E.S(a)  show  that  for  at 
least  the  first  five  cycles  of  the  flow  the  vortex  paths  are  virtually  identical.  Also  noted 
in  this  figure  are  sharp  cusps  where  the  overall  direction  of  the  vortices  is  abruptly 
chanaed. 


Although  the  difference  in  the  begining  locations  of  the  sixth  cycle  is  virtually 
inperceptible  (£  =  (0.995,  0.723)  vice  £  =  (0.970,  0.753)},  the  vortex  in  Figure  E.7(b) 
forms  a  sharp  cusp  at  the  top  of  the  cylinder  (geographical  location  of  the  maximum 
velocity  at  any  time  in  the  flow)  but  the  vortex  in  Figure  E.S(b)  is  captured  by  its 
image  and  orbits  the  cylinder.  The  ending  points  for  cycle  six  are  grossly  different  (see 
Figures  E.7(c)  and  E.S(c)}  and  the  vortex  paths  never  again  become  similar.  The  finite 
perturbation  in  starting  location  has  yielded  explicitly  different  results.  Nevertheless, 
the  vortex  paths  in  Figures  E.7(c)  and  E.S(c)  remain  bounded.  Extending  the 
calculations  to  100  cycles  of  the  flow.  Figure  E.9  shows  that  while  neither  vortex  paths 
are  identical,  both  remain  bounded.  They  clearly  seem  to  slide  in  and  out  of  chaotic 
behavior,  preceded  and  followed  by  periods  of  stability. 

Figure  E.10  is  a  plot  of  the  magnitude  of  the  velocity  of  the  vortices.  It  is  again 
evident  that  these  two  vortices  have  virtually  identical  velocity  magnitudes  until  the 
sixth  cycle  of  flow.  The  capture  and  looping  depicted  in  Figures  E.Tb)  and  E . S( b >  is 
also  evident  in  Figure  E.10.  Bieaking  down  this  velocity  magnitude  into  its  u  and  v 
components.  Figures  E.ll  through  E.13  shew  that  overall  periods  of  stability  are 
intermixed  with  these  vortex  captures. 

Figures  E.14  through  E.19  and  E.20  through  E.25  are  similar  investigations  that 
document  two  starting  points  of  the  vortices  where  both  initial  x-locations  were  set  at 
x  =  0.0  and  the  oscillating  flow  velocity  was  set  at  L:0  =  U+ax  and  U0  =  U“ax, 
respectively,  Figures  E.26  through  E.31  are  similar  investigations  that  document  two 
starting  points  of  the  vortex  in  the  vicinity  of  x=  1.0,  y=  1.0,  and  the  oscillating  flow 
velocity  was  set  at  L'o  =  0+.  These  investigations  support  the  conclusions  of  the  first 
investigation  and  demonstrate  once  again  a  finite  precission  input  does  not  yield  a 
finite  precission  output.  They  also  demonstrate  the  periods  of  stability  that  precede 
and  follow  chaotic  behavior,  (i.e.,  chaos  in  order  and  order  in  chaos). 

As  seen  in  the  proceeding  set  of  experiments,  the  looping  or  continued  looping  of 
the  vortex  tightly  around  the  cylinder,  was  a  function  not  only  of  the  location  of  the 
vortex  but  also  of  the  current  state  of  the  oscillating  flow.  This  suggests  that  at  a 
given  initial  location,  a  variance  of  initial  velocity  or  phase  shift  of  the  oscillating  flow 
could  drive  the  system  from  stable  to  chaotic  behavior  at  a  given  point  in  time.  Figure 
E.32  is  a  series  of  vortex  paths  that  supports  this  concept.  Four  vortex  paths  were 
generated  with  the  initial  location  of  x=y=  1.0  and  a  variance  of  initial  velocity  by  a 
phase  shift  from  <>.~5  |l  =  l'~  . )  to  1.0  (l.  =  iV1~  ).  It  is  clear  that  while  the  Fisurc 
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E. 32(a)  documents  stability  even  after  ten  cycles,  the  cusps  and  asymmetric  tracks  in 
Figure  E. 32(d)  demonstrate  the  mixture  of  stable  and  chaotic  regions. 

Figure  E.33  depicts  the  results  of  a  similar  investigation  where  the  initial  location 
of  the  vortices  is  at  (0.  2.2192)  and  (0,  2.2193).  These  results  also  show  that  the 
chaotic  behavior  depends  not  only  on  the  vortex  location  but  also  on  the  Dow 
environment.  In  fact,  the  two  are  coupled  inseparably. 

The  stability  and  symmetry  of  Figures  E. 32(a)  and  E. 33(a)  also  suggest,  however, 
that  chaotic  motion  may  not  be  inevitable.  They  infer  that  the  location  of  the  vortex 
and  field  How  may  be  tuned  in  such  a  way  as  to  permanently  avoid  chaos.  This  would 
be  paramount  to  permanently  balancing  on  a  razor's  edge.  However,  such  a  balance  is 
unstable.  As  seen  in  the  first  investigation,  the  slightest  disturbance  is  sufficient  to 
push  the  system  into  chaotic  behavior. 

The  final  phase  of  the  one  vortex  experiment  was  to  examine  the  effects  of  higher 
order  models  of  a  potential  vortex  in  an  effort  to  better  understand  vortices  in  nature. 
Figure  E.34  is  a  comparison  of  the  results  obtained  with  ideal,  Rankine,  and  Gaussian 
time  dependent  vortex  models.  It  is  clear  from  these  figures  that  while  the  onset  of 
chaos  may  be  different,  the  chaotic  behavior  and  susceptibility  to  chaos  remains 
virtually  the  same. 

C.  TWO  VORTEX  MODEL 

Separated  flow  past  a  bluff  body  generates  pairs  of  counterrotating  vortices.  In 
both  steady  and  oscillating  flow  about  a  circular  cylinder,  the  first  vortex  pair  is  formed 
symmetrically  and  is  significantly  stronger  than  subsequent  asymmetric  vortices. 
Additionally,  it  has  been  documented  that  in  oscillating  flow  in  the  region  S  <  K  <  13. 
a  half  von-Karman  vortex  street  forms  perpendicular  or  transverse  to  the  flow 
direction.  In  this  section  of  the  experimentation,  two  vortices  were  placed  in  the 
vicinity  of  a  circular  cylinder  immersed  in  a  harmonically  oscillating  flow.  The  vortices 
were  tracked  through  various  numbers  of  cycles  of  flow.  The  locations  and  strengths 
of  the  vortices  were  varied  in  order  to  examine  the  mechanism  for  the  formation  of  a 
half  von-Karman  street.  Additionally,  the  stability  of  the  vortices  as  discussed  in  the 
one  vortex  model  was  examined.  For  all  of  the  investigations,  the  characteristic  of  the 
flow  was  set  at  K  =10.0.  The  strength  of  the  vortices  was  set  at  K  =  0.50±6 

c  - 

(although  contrary  in  direction),  and  their  locations  at  C^  =  (a.a)  and 

=  i  a  ±  6,-a  ±  0  i.  where  O.SO<  a<  1.414,  and  6  <  <  1.  Based  on  flow  visualization,  the 

initial  velocitv  was  set  at  U  =  L ”... 

o  max 


The  first  part  of  the  investigation  required  the  analysis  of  vortices,  symmetric  in 
strength  and  location.  Figure  E.35  demonstrates  that  within  the  regions  examined,  the 
mutual  induction  velocities  between  the  real  vortices  dominate  the  advection  of  the 
vortices.  It  is  also  clear  that  the  vortices  are  quickly  convected  away  in  parallel  paths. 
The  initial  radial  location  determines  the  separation  distance  of  the  two  parallel  paths. 
This  separation  distance  combined  with  the  strength  of  the  vortices  determines  the  rate 
at  which  the  vortices  are  swept  away. 

Figures  E.36  through  E.3S  are  of  a  senes  of  investigations  where  the  vortices 
have  equal  strengths  (K^=Kg)  but  the  locations  are  asymmetric  f~[f  =  crx  2  a  i- 
The  overall  result  is  a  general  curving  of  the  vortex  path  toward  the  vortex  that  had 
the  larger  initial  radial  distance  from  the  cylinder.  This  is  a  result  of  the  image  vortices 
having  a  greater  effect  on  the  vortex  that  was  closer  to  the  cylinder.  It  is  seen  from 
Figure  E.3"’  that  a  10%  difference  (c  =0.90)  in  the  two  starting  locations  is  sufficient 
to  significantly  curve  the  vortex  paths.  The  vortex  paths  do  not  gravitate  to  each  other 
though  and  the  vortices  do  not  become  coincident. 

Figures  E.39  through  E.42  are  of  a  series  of  investigations  where  the  vortices 
have  symmetric  initial  locations  =  and  asymmetric  strengths  (Kg  =  cK  *  K^). 
The  result  here  is  a  curvature  similar  to  that  of  the  previous  investigation.  From 
Figure  E .4 1  it  is  seen  that  only  a  5%  difference  ( cK  =  0.95)  in  the  vortex  strengths  is 
sufficient  to  significantly  alter  the  vortex  paths.  Clearly  the  results  indicate  that  the 
system  is  much  more  sensitive  to  asymmetry  in  vortex  strengths  than  to  asymmetry  in 
vortex  location.  Again  the  vortices  do  not  become  coincident. 

Figure  E .43  is  a  combination  of  the  previous  investigations  where  there  is 
asymmetry  in  both  the  vortex  strengths  (Kg=cK*K^)  and  in  initial  location 
( |qg!  =  £r  *  |C ).  Again  both  vortex  tracks  bend  toward  the  dominant  vortex  but  as 
expected,  with  a  more  severe  curvature.  Although  in  real  flow  two  vortices  are 
generated  even.-  time  the  oscillating  flow  direction  changes,  this  studv  of  two  vortices 
moving  together  to  either  the  top  or  bottom  of  the  cylinder  is  particularly  significant. 
It  suggests  the  evolution  of  a  half  von-Karman  vortex  street  transverse  to  the  flow 
based  on  the  asymmetry  of  the  vortices  in  the  street. 

1  inallv,  an  investigation  was  made  where  the  asymmetry  of  the  vortices  was  so 
extensive  that  the  vortex  paths  curved  back  to  the  vicinity  of  the  cylinder.  Although  in 
viscous  flow  the  life  of  a  vortex  is  finite,  combinations  of  additional  vortices  coupled 
with  the  flow  field  make  this  situation  plausible  and  warrants  its  investigation.  Figure 


E.44  documents  the  sharp  cusps,  capture,  and  chaotic  behavior  of  the  vortices  as 
discussed  in  the  previous  section.  With  the  additional  vortex  (and  image)  the  coupling 
is  much  more  complicated  between  the  geographical  location  of  the  vertex  and  the 
flow  field.  The  identification  of  chaotic  behavior  is  none  the  less  possible. 


V.  CONCLUSIONS 


The  results  presented  herein  warranted  the  following  conclusions: 

1.  Finite  precision  information  about  the  initial  characteristics  of  a  vortex  (or  of 
larger  number  of  vortices)  provides  no  finite-precision  information  about  the 
characteristics  of  flow  at  later  times  and  that  the  lack  of  spanwise  coherence  of 
the  vortices  is  not  the  sole  cause  of  cycle-to-cycle  variations  in  the  measured 


4. 


6. 


quantities. 

An  idealized  vortex  model  suggests  that  the  sensitivity  of  the  characteristics  of 
the  flow  to  random  disturbances  superimposed  on  the  motion  of  vortices  at  an 
earlier  time  is  the  main  flow  feature  responsible  for  the  observed  cycle-to-cycle 
variations. 

The  chaotic  behavior  of  the  vortex  in  oscillating  flow  occurs  in  regions  of 
susceptibility  similar  to  those  found  by  Dowell  and  Pezeshki  [Ref.  4|  with 
similar  sensitivity . 

The  susceptible  regions  for  chaos  in  oscillating  flow  are  coupled  inseparably 
with  the  location  of  the  vortices  and  the  state  of  the  environmental  flow  field. 
In  other  worlds,  for  a  g;\ en  ;low  phase  there  is  a  vortex  location  that  will  lead 
to  chaos.  Likewise,  for  every  vertex  position  near  the  cylinder,  there  is  a 
particular  How  phase  that  w-il  ev : ntually  lead  to  chaos. 

The  formation  of  a  half  con- Karma::  vortex  street  in  the  ranee  of  8  <  K  <13 

w  V. 

may  be  attributed  to  the  uxymmet.-.  m  the  formation  of  the  vortices. 

The  extent  to  which  the  sensitivirv  of  the  flow  to  the  conditions  prevailing  in 
the  previous  cycles  is  damped  by  finite  core  effects  and  viscosity  has  been 


investieated  and  it  has  been  shown  that  thev  chance  or.lv  the  time  of  onset  of 


the  chaos,  not  its  mv evitable  occurrence 
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APPENDIX  A 

ONE  VORTEX  COMPUTER  MODEL 


C  *  * 

C  *  VORTEX  * 

C  *  * 

C  *  BY  LT.  W.T.  MCCOY  * 

C  *  * 

Q  ■k'k'k-k-k-k-k-k^'k-k-k'k'k-k'K-k-k^-k-k-^'k-k-k'k-k-k^-k^rk'k'k'kir-k-k-k'k-k-k'k^^-k-k-k-k^'A-k'k'k^^-k^'k'k-k'k^^rk-k-k-k'k 

c 

C  DESCRIPTION 

C 

C  THIS  PROGRAM  TRACKS  A  VORTEX  IN  COMPLEX  SPACE  AROUND  A  UNIT  CYLINDER 

C  IN  OSCILLATING  FLOW.  INITIAL  LOCATION,  STARTING  TIME,  AND  THE  NUMBER 

C  OF  CYCLES  OF  THE  OSCILLATING  FLOW  ARE  REQUIRED  INPUT. 

C 

C  THE  PROGRAM  IS  DESIGNED  TO  RUN  ON  DISSPLA  BY  USING  THE  FOLLOWING 
C  DISSPLA  EXEC  COMMAND : 1 DISSPLA  VORTEX’.  IT  IS  CURRENTLY 
C  DIMENSIONED  TO  RUN  A  MAX  OF  10  CYCLES  AT  1440  STEPS  PER  CYCLE. 

C  2.0M  OF  VIRTUAL  MEMORY  IS  REQUIRED  TO  RUN  THE  PROGRAM  AT  MAXIMUM 
C  CAPACITY.  THE  PROGRAM  IS  INTERACTIVE  FOR  INPUT  OPTIONS. 

C 

C  ALL  CONSTANTS  ARE  SET  IN  THE  MAIN  PROGRAM.  PREVIOUS  INPUT  OPTIONS 
C  ARE  READ  FROM  FILE  'FILE  FT25F001'.  FAILURE  TO  GENERATE  THIS  FILE 
C  PRIOR  TO  INITIAL  RUNNING  RESULTS  IN  A  ERROR  MESSAGE  AND  TERMINATION 
C  OF  THE  RUN.  PRIOR  TO  RUNNING,  FILE  'FILE  FT25F0011  SHOULD  5E 
C  INITIATED  CONTAINING  4  COMPLEX  NUMBERS.  1  REAL  NUMBER,  AND  I  INTEGER. 
C  WITH  THIS  INPUT  THE  PROGRAM  DOES  NOT  NEED  TO  BE  TO  BE  RECOMPILED 
C  FOR  EACH  RUN. 

C 

C  THIS  PROGRAM  IS  WRITTEN  IN  DOUBLE  PRECISION  AND  EMPLOYS  THE  MILNE 
C  PREDICTOR-CORRECTOR  CONVECTION  SCHEME,  WITH  INITIAL  POINTS  GENER- 
C  ATED  BY  A  4TH  ORDER  RUNGE-KUTTA  SCHEME. 


c 


C  OUTPUT  OF  THE  PROGRAM  INCLUDES  THE  LAST  FOUR  VALUES  OF  THE  LOCATION 
C  OF  THE  VORTEX,  FINISH  TIME,  AND  THE  NUMBER  OF  CYCLES.  THESE  VALUES 
C  ARE  WRITTEN  IN  FILE  'FILE  FT25F001'  AND  ALLOW  FOR  OBSERVATIONS  BEYOND 
C  THE  MAX  LIMIT  OF  10  CYCLES  BY  CONTINUING  WITH  ADDITIONAL  RUNS. 

C 

C  OUTPUT  FOR  THE  VORTEX  PATH  AS  WELL  AS  FORCE  COEFFICIENTS  AND 
C  VELOCITY  HUN I PULAT IONS  ARE  HADE  GRAPHICALLY.  SELECTIONS  OF  PLOTS 
C  FOR  OUTPUT  ARE  MADE  INTERACTIVELY.  PLOTS  MAY  BE  MADE  ON 
C  A  TEK613  OR  SENT  TO  A  METAFILE  (PRE-EDIT  SUBROUTINE  'GRAPH'). 

C 

C 

C . 

c 

C  ALPHABETICAL  LISTING  OF  VARIABLES 

C 

C  VARIABLE  (SUBROUTINE  USED)  DESCRIPTION. 


C 

C . 

c 

C  CD(N) . (GRAPH) NON-DIMENSIONAL  ARRAY  OF  INLINE  (DRAG) 

C  FORCE  COEFFICIENTS  FOR  VORTEX  POSITION  Z(N) . 

C 

C  CL ( N ) . (GRAPH) NON -DIMENSIONAL  ARRAY  OF  PERPENDICTULAR 

C  (LIFT)  FORCE  COEFFICIENTS  FOR  VORTEX  POSITION 

C  Z (N) . 

C 

C  -  COEF . (FORCE, GRAPH) COMPLEX  VARIABLE  CONTAINING  AS  ITS 

C  PARTS  THE  NON-DIMENSIONAL  FORCE  COEFFICIENTS. 

v- 

C  DELTAT . ( STEP, Q, OUTPUT, FORCE, GRAPH) INCREMENTAL  STEP 

C  IN  NON-DIMENSIONAL  TIME  (INVERSE  OF  STEPS  PER 

C  CYCLE ) . 

C 

C  DELZ . (STEP) COMPLEX  INCREMENTAL  STEP  OF  THE  LOCATION 


V 
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OF  THE  VORTEX  (FROM  Z1  TO  Z2) . 


DFDZ . (Q) COMPLEX  PARTIAL  DERIVATIVE  OF  THE  COMPLEX 

FUNCTION  F  WITH  TO  RESPECT  Z. 

DISP(N) . (STEP, GRAPH) ARRAY  CONTAINING  THE  TOTAL  DISPLACEMENT 

OF  THE  VORTEX  FOR  A  CORRESPONDING  LOCATION  Z(N). 

I . (Q, FORCE) CONSTANT  COMPLEX  NUMBER  =  (0.0, 1.0) 

ILHS . (GRAPH) LEFT  HAND  VALUE  OF  THE  TIME  AXIS  OF  THE 

FORCE  COEFFICIENT  PLOT. 

IRHS . (GRAPH) RIGHT  HAND  VALUE  OF  THE  TIME  AXIS  OF  THE 

FORCE  COEFFICIENT  PLOT. 

ISTEPS . (ALL) NUMBER  OF  INCREMENTAL  STEPS  PER  CYCLE  TO  BE 

TAKEN  FOR  ALL  ITERATIONS. 

Kl-4 . (RUNGE)  RUNGE-KUTTA  CONSTANTS 

KAPPA . (ALL) STRENGTH  OF  THE  VORTEX  (GAMMA/ (2*?I*UMAX) ) 

KC . (ALL)  KEULLEGAN- CARPENTER  NUMBER  FOR  THE  OSCILLATING 

FLOW. 

LI . (INPUT, GRAPH)  LOGICAL  OPERATOR. 

MUMCYC . ( INPUT) TOTAL  NUMBER  OF  CYCLES  TO  BE  EVALUATED. 

NUMPTS . (ALL)TOTAL  NUMBER  OF  POINTS  TO  BE  DETERMINED. 

(MUMCYC  *  ISTEPS) 

PHSHFT . (GRPAH) INITIAL  PHASE  SHIFT  OF  THE  FLOW  (0<PHSHFT<1) 


IWI*1HI*'*!L* 


C  PI..- . (ALL)  CONSTANT  =  3.141592654 

C 

C  Q . (Q) FUNCTION  SUBROUTINE  WHERE  Q  IS  THE  COMPLEX 

C  VELOCITY  OF  THE  VORTEX. 

C 

C  QMAG(N) . (GRAPH) ARRAY  CONTAINING  THE  MAGNITUDE  OF  THE 

C  COMPLEX  VELOCITY  QVECT(N) . 

C 

C  QVECT(N) . (STEP, GRAPH) COMPLEX  ARRAY  CONTAINING  THE  VELOCITY 

C  OF  THE  VORTEX  AT  ANY  POINT  Z(N). 

C 

C  RV . (Q)  VECTOR  BETWEEN  THE  VORTEX  AND  ITS  IMAGE 

C  WITHIN  THE  CYLINDER 

C 

C  T . ( GRAPH )  NON  -  DIMENSIONAL  TIME  (0-1  FOR  A  CYCLE). 

C 

C  TE . (RUNGE  ,  Q)NON-DIMENSIONAL  ENVIRONMENTAL  FLOW  TIME. 

C 

C  TT . ( STEP, Q,  FORCE) NON-DIMENSIONAL  TIME  CHANGED  TO 

C  RADIANS  FOR  COS/SIN  EVALUATIONS  (0-2*PI  FOR  A  CYCLE). 

C 

C  TTF . (OUTPUT)NON-DIMENSIONAL  FINISH  TIME  OF  THE  RUNS. 

C 

C  TTI . (ALL) INITIAL  NON-DIMENSIONAL  STARTING  TIME  FOR  THE 

C  OSCILLATING  FLOW. 

C 

C  TV . (Q)NON- DIMENSIONAL  TIME  LIFE  OF  THE  VORTEX. 

C  (TE  +  CONSTANT  FOR  R-CORE  AT  INCEPTION) 

C 

C  U . ( GRAPH  )X- DIRECT  I  ON  VELOCITY  MAGNITUDE. 

C 

C  V . ( GRAPH )Y- DIRECTION  VELOCITY  MAGNITUDE. 

C 

C  Z . (ALL)  INSTANTANEOUS  COMPLEX  LOCATION  OF  THE  VORTEX. 

C 


30 


C  ZI . (GRAPH)  IMAGINARY (Y)  LOCATIONOF  THE  VORTEX. 

C 

C  ZP . (Q) COMPLEX  PREDICTED  VALUE  FOR  Z(N) 

C 

C  ZR . (GRAPH) REAL (X)  LOCATION  OF  THE  VORTEX. 

C 

C  ZZ . (Q) DUMMY  ARGUMENT  FOR  THE  FUNCTION  SUBROUTINE  Q 

C  WHICH  CONTAINS  THE  PASSED  VALUE  OF  Z(N). 

C  (RUNGE)  DUMMY  ARGUMENT  CALL  FUNCTION  SUBROUITME  Q 

C 

C  . 

c 

c 


(2 

C  *  MAIN  PROGRAM  * 

c 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  KC, KAPPA 

REAL*4  DISP ( 14409 ) , CD( 14409 ), CL( 14409 ) 

COMPLEX* 16  Z ( 14409 ) ,QVECT( 14409) 

COMMON  /FIRST/  TTI , ISTEP5  /SECOND/  NUMPTS  /THIRD/  KC, KAPPA 
1  /FOURTH/  PI 

C 

C  SET  CONSTANTS 
C 

KAPPA  =  0.5D0 
KC  =  10. DO 
PI  =  3 . 141592654D0 
ISTEPS  =  1440 
C 

10  CALL  INPUT  (Z) 

CALL  STEP  ( Z , QVECT , DI SP ) 

CALL  OUTPUT (Z) 

CALL  FORCE (Z , CD , CL ) 


CALL  GRAPH  (Z , QVECT , DISP , CD , CL) 

C 

WRITE  (6,*) 'DO  YOU  WISH  TO  RE-RUN  THE  PROGRAM?  1-YES,  O-NO1 
READ  (5,*)  LI 
IF  (Ll.NE.l)  GOTO  20 

WRITE  (&,*)  'EXECUTION  CONTINUES  .  ' 

GOTO  10 
20  CONTINUE 
C 

CALL  DONEPL 
STOP 
END 
C 
C 

Q  *****************************************************  **********^*'k**‘^ 

C  *  SUBROUTINE  INPUT  * 

( 2  ********************************************** *****************'X***** 

c 

C  THIS  SUBROUTINE  INTERACTIVELY  INPUTS  THE  VALUES  FOR  THE  INITIAL 
C  LOCATION  OF  THE  VORTEX,  STARTING  TIME,  AND  THE  NUMBER  OF  CYCLES 
C  TO  BE  EXAMINED. 

C 

SUBROUTINE  INPUT  (Z) 

IMPLICIT  REALMS  (A-H,0-Z) 

COMPLEX* 16  Z( 14409) 

COMMON  /FIRST/  TTI , ISTEPS  /SECOND/  NUMPTS 
C 

C  READ  THE  LAST  4  POINTS  FOR  CONTINUATION  OPTION 
C 

REWIND  25 
DO  10  M  =  1,4 
10  READ (25 , *)  Z (N ) 

READ (25,*)  TTI 
READ ( 25 , * )  NUN CYC 
C 


C  CONTINUATION  OPTION 
C 

WRITE (6,*)  ’DO  YOU  WISH  TO  CONTINUE  THIS  RUN  FORM  THE  LAST?  1- 
#  0-NO1 

READ (5 , *)  LI 
IF  (L1.NE.0)  GOTO  20 

WRITE ( 6, *) 1  INPUT  THE  NEW  STARTING  POINT  Z(l):  (REAL , IMAG) 1 
READ (5,*)  Z(l) 

WRITE (6,*) 'INPUT  THE  NEW  VALUE  OF  TTI : 1 
READ ( 5 , * )  TTI 
C 

C  CALL  SUBROUTINE  RUNGE  TO  GENERATE  THE  FIRST  4  POINTS 
C 

CALL  RUNGE  (Z) 

C 

C  INPUT  THE  NUMBER  OF  CYCLES  TO  BE  ANALYZED 
C 

20  WRITE (6,*) 'THE  LAST  VALUE  FOR  NUMCYC  =', NUMCYC 

WRITE (6 , *) 'DO  YOU  WISH  TO  CHANGE  IT?  1-YES  O-NO1 
READ (5,*)  LI 
IF (LI. ME. 1)  GOTO  30 

WRITE (6,*)  'INPUT  THE  NEW  VALUE  OF  NUMCYC:' 

READ (5,*)  NUMCYC 
30  CONTINUE 
C 

NUMPTS  =  NUMCYC  *  ISTEPS  +  4 
C 

RETURN 

END 

C 

r 

Q  ■k'k-K7r-k-k~-i'7'-h-h'k-r:'k-k-k+:-k-k-ki'-K-r’k-'r-H-k-X'n*:  +  +j*-k-k7'.Tr-h-ry'-ki<-k-k-k-k-k-k-k-k-k-k-k-k'k-k’k-k-k-k'k7<'k'k  +  -k'kr: 

C  *  SUBROUTINE  RUNGE 


C  THIS  SUBROUTINE  USES  A  4TH  ORDER  RUNGE-KUTTA  METHOD  TO  GENERATE  THE 
C  FIRST  FOUR  VALUES  OF  Z. 

C 

SUBROUTINE  RUNGE  (Z) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEX* 16  Z( 14409) ,Q,K1 , K2 , K3 , K4 , ZP1 , ZZ 
COMMON  /FIRST/  TTI,ISTEPS 
C 

DELTAT  =  1 .DO/DFLOAT(ISTEPS) 

C 

C  RUNGE-KUTTA  EVALUATION  FOR  FIRST  FOUR  POINTS 
C 

DO  5  N  =  1,3 
ZZ  =  Z(N) 

TE1  =  (DFLOAT(N-l)  *  DELTAT  +  TTI ) 

K1  =  DLETAT  *  Q(ZZ,TE1) 

ZZ  =  Z(N)+0.5D0*K1 
TE2  =  TE1+  0 . SDO*DELTAT 
K2  =  DELTAT  *  Q(ZZ,TE2) 

ZZ  =  Z(N)+0 . 5DO*K2 

K3  =  DELTAT  *  Q(ZZ,  TE2 ) 

ZZ  =  Z(N)  +  K3 

TE4  =  TE1  +  DELTAT 

K4  =  DELTAT  *  Q(ZZ,  TE4 ) 

C 

5  Z(N+1)  =  Z(N)  +  (Kl-2 . D0*K2+2 . D0*K3+K4 ) /6 . DO 
C 
C 
C 

RETURN 

END 

C 

c 

C  *  SUBROUTINE  STEP  * 
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Q  ■k'k-k-k*-k'k,k'k'fi:'k'k'k’k*-rrk-k  -k-k  ********** 

c 

C  THIS  SUBROUTINE  PROGRESSIVELY  CALCULATES  THE  VALUES  OF  Z(N),  BASED  ON 
C  A  PREDICTOR-CORRECTOR  CONVECTION  SCHEME. 

C  IT  ALSO  DETERMINES  THE  VELOCITY  VECTOR,  MAGNITUDE  AND  PARTS  FOR 
C  FUTURE  MUNIPULATION  AS  WELL  AS  OVERALL  DISPLACEMENT  MAGNITUDES. 

C 

SUBROUTINE  STEP  (Z , QVECT , DISF ) 

IMPLICIT  REALMS  (A-H,0-Z) 

REAL*4  DISP (14409 ) 

REAL*3  KC, KAPPA 

COMPLEXES  Z(  14409)  ,  QVECT  ( 14409 )  ,  DELZ ,  ZP  ,  Q  ,  ZZ 

COMMON  /FIRST/  TTI,ISTEPS  /SECOND/  NUMPTS/THIRD/  KC, KAPPA 

r 

DELTAT  =  1 .DO/DFLOAT(ISTEPS) 

C 

C  EVALUATE  VELOCITIES  AND  DISPLACEMENTS  FOR  FIRST  4  POINTS 
C 

DISP ( 1 )  =  0.0 
C 

DO  7  N  =  1,4 

TE  =  (DFLOAT(N-l)  *  DELTAT  +  TTI ) 

ZZ  =  Z ( N ) 

QVECT (M)  =  Q ( ZZ , TE ) 

7  DISP (N+l )  =  Z(N+1)  -  Z(M) 

C 

DO  10  N  =  5,  NUMPTS 

NM1  =  N-l 

MM2  =  N-2 

NM3  =  N-3 

NM4  =  N-4 

TE  =  (DFLOAT (N-l ) *DELTAT+TTI ) 

C 

C  CALCULATE  PREDICTED  VALUE  OF  Z(N) 


ZP  =  Z(NM1)+(55.D0*QVECT(NM1)  -59 .D0*QVECT(NM2)  +37 . DO*QVECT (NM3 ) 
1  -9 . D0*QVECT(NM4) )*DELTAT*  2. DO  *  KC  /  24. DO 


\-V  >  V 


CALCULATE  ACTUAL  VALUE  OF  Z(N) 


DELZ  =  ( 9 . D0*Q (ZP , TE )  +  19 . D0*QVECT (NM1 )  -5 . D0*QVECT (NM2 ) 
1  +OVECT (MM3 )  ) *DELTAT  *  2. DO  *  KC  /  24. DO 


Z(N)  =  Z(NK1)  +  DELZ 


EVALUATE  VELOCITIES  AT  Z(N) 


QVECT(N)  =  Q(Z(N),TE) 

10  DISP(N)  =  DISP(NMl)  +  ABS (DELZ) 


RETURN 


FUNCTION  SUBPROGRAM  Q 


THIS  FUNCTION  SUBPORGRAM  IS  USED  TO  EVALUATE  VALUES  OF  Q 


COMPLEX  FUNCTION  Q*16(ZZ,TE) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*8  KC,  KAPPA 
C0MPLEX*16  ZZ,  DFDZ ,  I , RV 

COMMON  /FIRST/  TTI,ISTEPS  /THIRD/  KC, KAPPA  /FOURTH/  PI 


I  =  (0 . DO , 1 . DO  ) 


DELTAT  =  1  .DO/DFLOAT(ISTEPS) 


TT  =  TE  *  2. DO*  PI 


RV  =  2Z  -  1 . DO  /  DCONJG(ZZ) 


C 

DFDZ  =DSIN(TT)*(1 .DO-l.DO/(ZZ*ZZ) )-I*KAPPA/RV 
C 

Q  =  DCONJG(DFDZ) 

C 

RETURN 

END 

C 

c 

C  **************x********x*x*x*x*J<*****x***'k'k,k‘k'k'k-k-k*'k'k*:-kx-k'k-k-k-k-k-k-k-k-k'k-k-k-k 

C  *  SUBROUTINE  OUTPUT  * 

Q  ^^^’k’k'k'k'k'k'k'ky^'k-k’k'k-k'k-k-k'k-k-k-k-k'k-^-k-k’k'k'kTr-k’k-k-k'k’k-k-k^'k’k'k'k-k-ky^'k'k-k-kir'k'k'k'k'k'k-k-k'k'k'kTr’k'k'k 

c 

C  THIS  SUBROUTINE  IS  USED  TO  SAVE  THE  FINAL  VALUES  OF  Z,TT,  AND  NUMCYC 
C  FOR  POSSIBLE  CONTINUATION  OF  THE  PLOTS. 

C  THE  VALUES  ARE  SAVED  IN  FILE  ‘FILE  FT25F001'. 

C 

SUBROUTINE  OUTPUT (Z) 

IMPLICIT  REAL*8  (A-H,  O-Z) 

COMPLEX*16  Z( 14409) 

COMMON  /FIRST/  TTI,ISTEPS  /SECOND/  NUMPTS 
C 

DELTAT  =  1 . DO /FLOAT ( I STEPS ) 

NN  =  NUMPTS  -  4 
TTF  =  DFLOAT (NN)*  DELTAT  +  TTI 
NUMCYC  =  NN  /  I STEPS 
C 

REWIND  25 
DO  10  N  =  1 ,4 
NN  =  NN  +  1 
10  WRITE (25 , *)  Z ( NN ) 


WRITE (2 5,*)  TTF 
WRITE (2 5,*)  NUMCYC 


RETURN 


END 

C 

C 

C  *  SUBROUTINE  FORCE 

c 

C  THIS  SUBROUTINE  IS  USED  TO  CALCULATE  THE  FORCE  COEFFICIENTS  CD,  AND 
C  CL  FOR  THE  CYLINDER  BASED  ON  THE  LOCATIONS  OF  THE  VORTEX  AND  A 
C  FINITE  DIFFERENCING  SCHEME. 

C 

SUBROUTINE  FORCE (2 , CD , CL ) 

IMPLICIT  REAL'S  (A-H,  C-Z) 

REAL*4  CD ( 1440 9 ) , CL( 14409 ) 

REALMS  KC,  KAPPA 
COMPLEX'S  COEF 
COMPLEXES  Z(  14409),  I 

COMMON  /FIRST/  TTI.ISTEPS  /SECOND/  NUMPTS  /THIRD/  KC, KAPPA 
1  /FOURTH/  PI 

C 

I  =  ( 0  .  DO ,  1 . DO ) 

DELTAT  =  1  .DO/DFLOAT(ISTEPS) 

C 

DO  10  N  =  2,  NUMPTS 
C 

TT  =  (DFLOAT(N-l)  *  DELTAT  +  TTI)  *  2. DO  *  PI 


COEF  =  ( 2 . D0'PI*?I/KC) ^DCOS (TT) 

1  +  ( (PI'KAPPA)/ (KC'DELTAT) )*I* 

2  ( (Z(N)-l .DO/DCCNJG(Z(M) ) )- 

3  (Z(N-l)  -  1 . DO/ DCOM JG (Z(N-l)))) 
COEF  =  CONJG ( COEF ) 

CD(N)  =  REAL ( COEF ) 
iO  CL ( N )  =  I MAG ( COEF ) 


3S 


RETURN 


END 


C 

r 


SUBROUTINE  GRAPH 


SUBROUTINE  GRAPH  (Z , QVECT , DIS? , CD , CL) 

IMPLICIT  REALMS  (A-H,  O-Z) 

REAL*4  ZR ( 14409 ) ,ZI(14409) , X( 361 ) , Y( 361 ) , QMAG( 14409 ) 
REAL *4  DISP ( 14409 ) , CD ( 14409 ) , CL ( 14409 ) , T ( 14409 ) 

REAL~'4  U(  14409)  ,v  (14409)  , GARRAY ( 1 ) 

COMPLEX* 16  Z( 14409) , QVECT ( 14409 ) 

COMMON  /FIRST/  TTI , ISTEPS  /SECOND/  NUMPTS  /FOURTH/  PI 
C 

C  DEFINE  CYLINDER  OF  RADIUS  1.0 
C 

DO  10  N  =  1,361 

RAD  =  FLOAT ( N- 1 )  *  PI  /I 30.0 

X(N)  =  COS (RAD) 

10  Y(N)  =  SIN (RAD) 

C 

C  DEFINE  REAL  SINGLE  PRECISION  NUMBERS  FOR  PLOT  ROUTINES 
C 

DELTA!  =  1 .DO/DFLOAT(ISTEPS) 

DO  20  N  =1,  NUMPTS 
ZR(N)  =CREAL(Z(N) ) 

ZI(N)  =DIMAG ( Z (N) ) 

QMAG(N)  =  ABS ( QVECT (M ) ) 

U(N)  =  DREAL (QVECT (N) ) 

20  V ( N )  =  D I HAG ( QVECT ( M ) ) 
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.  -V.N  -V. 


-■%•••/  V 


c 


C  GENERATE  PLOTS 
C 

CALL  COMPRS 
C  CALL  TEK613 

CALL  BLOWUP ( .3365655) 

C  CALL  BLOWUP (.55) 

C  CALL  BLOWUP (2.0) 

C  CALL  NOERDR 

C 

C  VORTEX  TRACK  PLOT 
C 

CALL  GRACE  (0.0) 

CALL  NOCHEK 

CALL  PAGE  (11.0,3.5) 

CALL  AREA2D(10. ,  5.71425) 

CALL  YNAME ( ' Y $  1 ,100) 

CALL  XNAME ( >  X$ 1 ,100) 

CALL  YAXANG (0.0) 

CALL  1NTAXS 
C 

CALL  XTICKS ( 2 ) 

CALL  YTICKS ( 2 ) 

CALL  GRAF (-7, 1,7, -4, 1,4) 

C 

CALL  CURVE(X,Y, 361,0) 

GARRAY ( 1 )  =0.03 

CALL  SHADE ( X , Y , 36 1 , 4 5 , GARRAY ,1,0,0) 

C 

CALL  CURVE (ZR,ZI ,NUMPTS , 2000C) 

C  CALL  CURVE (ZR, 21, NUMPTS,0) 

CALL  EMDPL(O) 

C 

C  HAKE  CONTINUATION  SELECTIONS 
C 

WRITE ( 6 , * )  '  CONTINUATION  NEMO 


WRITE (6,*) 

WRITE (6 ,  *) 

WRITE (6,*) 'TERMINATION  OF  PROGRAM . O' 

WRITE (6  ,*)  1  RERUN  PROGRAM . 1' 

WRITE (6 ,  *)  '  CONTINUE  WITH  NORMAL  SEQUENCE  OF  PLOTS . 2' 


WRITE ( o  ,  * ) 

WRITE (6 , *)  '  PLEASE  ENTER  YOUR  SELECTION:  (0,1,2)' 

READ (5,*) LI 
IF(Ll.LT.l)GOTO  100 
IF (LI . EQ . 1 )GOTO  100 
CONTINUE 
C 

PHSHFT  =  TTI  -  DFLOAT( INT(TTI ) ) 

C 

DO  40  N  =  1 , NUMPTS 

40  T(N)  =  (DFLOAT(N-l)  *  DELTAT  +  TTI  -  PHSHFT) 

C 

WRITE (6,*) 'DO  YOU  WANT  A  FORCE  PLOT?  1-YES,  0-N0' 

READ (5,*)  LI 
IF  (Ll.NE.l)GOTO  50 
C 

C  FORCE  VS  TIME  PLOTS 
C 

WRITE (6,*)  'ENTER  VALUE  OF  THE  LEFT  HAND  SIDE  OF  THE  TIME  SCALE: 
READ (5,*)  ILHS 

WRITE (6,*)  ’ENTER  VALUE  OF  THE  RIGHT  HAND  SIDE  OF  THE  TIME  SCALE 
READ ( 5 , * )  IRHS 
C 

CALL  AREA2D (10.0,  7.5) 

CALL  YNAME ( ’ COEFS 1 ,100) 

CALL  XNAME ( 'T*$ ' ,100) 

CALL  YAXANG (0.0) 

CALL  INTAXS 
C 

CALL  GRAF ( ILHS , . 5 , IRHS , -4 , 1 . , 4 . 0 ) 


4! 


CALL  CURVE (T, CD, NUMPTS, 0) 
CALL  DASH 

CALL  CURVE ( T , CL , NUMPTS , 9 ) 
CALL  RESET  ( 1  DASH'  ) 

CALL  EMDPL(O) 


50  WRITE (6 , * ) 'DO  YOU  WANT  A  Q-MAG  VS  DISP  PLOT?  1-YES,  0-N0 
READ (5 ,  *)  LI 
IF  (Ll.NE.l)GOTO  60 
C 

C  MAGNITUDE  OF  Q  VS.  DISPLACEMENT  PLOT 


CALL  AREA2D ( 10 . 0 ,  7.5) 
CALL  YNAME ( 1 Q$ 1 , 100) 
CALL  XNAME ( 1  S3 ' ,100) 
CALL  YAXANG (0.0) 

CALL  INTAXS 


CALL  GRACE (2.) 

CALL  GRAF(0, 10,160, 0,1. ,10.0) 

CALL  CURVE ( DISP, QMAG, NUMPTS, 0) 

CALL  ENDPL(O) 

60  WRITE (6,*) 'DO  YOU  WANT  A  Q-MAG  VS  TIME  PLOT’  1-YES.  0-N0 
READ (5,*)  LI 
IF  (Ll.NE.l)GOTO  70 

MAGNITUDE  OF  Q  VS .  TIME  PLOT 


CALL  AREA2D (10.0,  7.5) 
CALL  YNAME ( ’ Q3  * , 100 ) 
CALL  XNAME ( 1 T*$ 1 , 100 ) 
CALL  YAXANG (0.0) 


CALL  INTAXS 


C 

CALL  GRAF(0, 1,10, 0,1. ,10.0) 

C 

CALL  CURVE ( T , QMAG , NUMPTS , 0 ) 

CALL  EMDPL(O) 

C 

70  WRITE (6,*) ’DO  YOU  WANT  A  U-VEL  VS  X-POS  PLOT?  1-YES,  0-N0 
READ (5,*)  LI 
IF  (Ll.NS.l)GOTO  SO 
C 

C  U-VELOCITY  POSITION  PLOTS 
C 

CALL  AREA2D (10.0,  7.5) 

CALL  YNAME( 'U-VELOCITYS ' , 100) 

CALL  XNAME ( ' X-POSITION$ ' , 100 ) 

CALL  YAXANG (0.0) 

CALL  INTAXS 
C 

CALL  GRAF(-7,1,7,-10,2. ,10.0) 

C 

CALL  CURVE (2R,U, NUMPTS, 0) 

CALL  ENDPL(O) 

C 

80  WRITE (6 , *) 1  DO  YOU  WANT  A  V-VEL  VS  Y-POS  PLOT?  1-YES,  0-N0 
READ ( 5 , *)  LI 
IF  (Ll.NE.l)GOTO  90 
C 

C  V-VELOCITY  POSITION  PLOTS 
C 

CALL  AREA2D (10.0,  7.5) 

CALL  YNAHE ( 1 V-VELOCITYS ' , 100) 

CALL  XMAME ( ' Y-POSITIONS 1 ,100) 

CALL  YAXANG (0.0) 

CALL  INTAXS 


CALL  GRAF (-4, . 5 ,4, -7 , 1 . ,7.0) 

C 

CALL  CURVE (ZI , V,NUMPTS , 0 ) 

CALL  ENDPL(O) 

C 

90  WRITE (6 , *) 1  DO  YOU  WANT  A  U-VEL  VS  V-VEL  PLOT?  1-YES,  0-N0 
READ (5 , *)  LI 
IF  ( LI .ME . 1 ) GOTO  100 
C 

C  U  VS  V  VELOCITY  PLOTS 
C 

CALL  AREA2D (7 . 5 ,  7.5) 

CALL  YNAME ( ' V-VELOCITYS 1 ,100) 

CALL  XNAME ( ' U-VELOCITYS ' , 100) 

CALL  YAXANG (0.0) 

CALL  INTAXS 
C 

CALL  GRAF (-10,2,10,-10,2. ,10.0) 

C 

CALL  CURVE (U , V ,NUMPTS , 0) 

CALL  EMDPL(O) 

C 

100  CONTINUE 
C 


RETURN 


APPENDIX  B 

RANK1NE  VORTEX  COMPUTER  MODEL 


c 

C  THIS  FUNCTION  SUBROUTINE  IS  WRITTEN  FOR  USE  WITH  PROGRAM  'VORTEX 
C  FORTRAN ' .  IT  CONTAINS  VELOCITY  CALCULATIONS  FOR  A  RANKINE 
C  VORTEX  MODEL.  IT  MAY  BE  INSERTED  IN  PLACE  OF  THE  IDEAL 
C  VORTEX  MODEL  (FUNCTION  SUBROUITNE  Q  IN  'VORTEX  FORTRAN1) 

C 

C 


C  |  FUNCTION  SUBPROGRAM  Q  (RANKINE  MODEL) 


C 

C  THIS  FUNCTION  SUBPORGRAM  IS  USED  TO  EVALUATE  VALUES  OF  Q 
C 

COMPLEX  FUNCTION  Q*16(ZZ,TE) 

IMPLICIT  REAL *8  (A-H,0-Z) 

REAL*8  KC,  KAPPA 

COMPLEXES  ZZ,  DFDZ ,  I,  RV,  VINDUC 

COMMON  /FIRST/  TTI,ISTEPS  /THIRD/  KC , KAPPA  /FOURTH/  PI 
C 

RCORE  =  0.1D0 
I  =  ( 0 . DO , 1 . DO ) 

DELTAT  =  1 .DO/DFLOAT(ISTEPS) 

TT  =  TE  *  2. DO  *  PI 
C 

RV  =  ZZ  -  1.D0  /  DCONJG(ZZ) 

CHECK  —  ABS(RV) 

r 

IF  (RCORE .GT. CHECK)  GOTO  10 
C 


VINDUC  =  - 1 .  DO  /  RV 
GOTO  20 


c 

10  VINDUC  =  -RV  /  RCORE**2 
C 

20  DFDZ  =DSIN(TT)*(1  .DO-1  .DO/ (ZZ*ZZ) )+KAPPA*I*VINDUC 
C 

Q  =  DCOHJG(DFDZ) 

C 

RETURN 

END 


APPENDIX  C 

GAUSSIAN  VORTEX  COMPUTER  MODEL 


c 

C  THIS  FUNCTION  SUBROUTINE  IS  WRITTEN  FOR  USE  WITH  PROGRAM  'VORTEX 
C  FORTRAN1.  IT  CONTAINS  VELOCITY  CALCULATIONS  FOR  A  GAUSSIAN  TIME 
C  DEPENDENT  VORTEX  MODEL.  IT  MAY  BE  INSERTED  IN  PLACE  OF  THE  IDEAL 
C  VORTEX  MODEL  (FUNCTION  SUBROUITNE  Q  IN  'VORTEX  FORTRAN') 

C 
C 

C . . - . . . - . - . 

C  |  FUNCTION  SUBPROGRAM  Q  (GAUSSIAN  MODEL) 

C . - . - . 

c 

C  THIS  FUNCTION  SUBPORGRAM  IS  USED  TO  EVALUATE  VALUES  OF  Q 
C 

COMPLEX  FUNCTION  Q*16(ZZ,TE) 

IMPLICIT  REAL*8  (A-H.O-Z) 

REAL*3  KC,  KAPPA 
COMPLEX* 16  ZZ,  DFDZ ,  I,RV 

COMMON  /FIRST/  TTI,ISTEPS  /THIRD/  KC,  KAPPA  /FOURTH/  PI 
C 

I  =  ( 0 . DO , 1 . DO ) 

DELTAT  =  1 .DO/DFLOAT(ISTEPS) 

TT  =  TE  *  2. DO*  PI 
TVINV  =  1 , DO  /(TE  +  . 497440635D0 ) 

C 

RV  =  ZZ  -  1.D0  /  DCONJG(ZZ) 

RVM/iG  =  ABS(RV) 

ARG  =  -62 . 5D0*TVINV*RVMAG**2 

IF  (ARG. LE. 50. DO)  GOTO  10 
ARG  =  50. DO 

r 

10  VMOD  =  1 . DO  -  DEXP(ARG) 


c 

DFDZ  =DSIN(TT)*(1.DO-1.DO/(ZZ*ZZ))-I*(KAPPA/RV)*VMOD 
C 

Q  =  DCONJG(DFDZ) 

C 

RETURN 

END 
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APPENDIX  D 

TWO  VORTEX  COMPUTER  MODEL 


X  * 

*  TWOVOR  * 

*  *f 

*  BY  LT.  W.T.  MCCOY  * 

k  k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

DESCRIPTION 

THIS  PROGRAM  TRACKS  TWO  VORTICIES  IN  COMPLEX  SPACE  AROUND  A  UNIT 
CYLINDER  IN  OSCILLATING  FLOW.  INITIAL  LOCATION,  STARTING  TIME,  AND 
THE  NUMBER  OF  CYCLES  OF  THE  OSCILLATING  FLOW  ARE  REQUIRED  INPUT. 

THE  PROGRAM  IS  DESIGNED  TO  RUN  ON  DISSPLA  BY  USING  THE  FOLLOWING 
DISSPLA  EXEC  COMMAND DISSPLA  TWOVOR'.  IT  IS  CURRENTLY 
DIMENSIONED  TO  RUN  A  MAX  OF  6  CYCLES  AT  1440  STEPS  PER  CYCLE. 

2.0M  OF  VIRTUAL  MEMORY  IS  REQUIRED  TO  RUN  THE  PROGRAM  AT  MAXIMUM 
CAPACITY.  THE  PROGRAM  IS  INTERACTIVE  FOR  INPUT  OPTIONS. 

ALL  CONSTANTS  ARE  SET  IN  THE  MAIN  PROGRAM.  PREVIOUS  INPUT  OPTIONS 
ARE  READ  FROM  FILE  'FILE  FT35F001 1 .  FAILURE  TO  GENERATE  THIS  FILE 
PRIOR  TO  INITIAL  RUNNING  RESULTS  IN  A  ERROR  MESSAGE  AND  TERMINATION 
OF  THE  RUN.  PRIOR  TO  RUNNING,  FILE  'FILE  FT25F001 '  SHOULD  BE 
INITIATED  CONTAINING  3  COMPLEX  NUMBERS,  i  REAL  NUMBER,  AND  1  INTEGER. 
WITH  THIS  INPUT  THE  PROGRAM  DOES  NOT  NEED  TO  BE  TO  BE  RECOMPILED 
FOR  EACH  RUN. 

THIS  PROGRAM  IS  WRITTEN  IN  DOUBLE  PRECISION  AND  EMPLOYS  THE  MILNE 
PREDICTOR-CORRECTOR  CONVECTION  SCHEME,  WITH  INITIAL  POINTS  GENER¬ 
ATED  BY  A  4TH  ORDER  RUMGE-KUTTA  SCHEME. 


c 

C  OUTPUT  OF  THE  PROGRAM  INCLUDES  THE  LAST  FOUR  VALUES  OF  THE  LOCATION 
C  OF  THE  VORTICES,  FINISH  TIME,  AND  THE  NUM3ER  OF  CYCLES.  THESE  VALUES 
C  ARE  WRITTEN  IN  FILE  'FILE  FT35F001'  AND  ALLOW  FOR  OBSERVATIONS  BEYOND 
C  THE  MAX  LIMIT  OF  6  CYCLES  5Y  CONTINUING  WITH  ADDITIONAL  RUNS. 

C 

C  OUTPUT  FOR  THE  VORTEX  PATHS  ARE  MADE  GRAPHICALLY.  PLOTS  MAY  BE  MADE 
C  CM  A  TEX518  OR  SENT  TO  A  METAFILE  (PRE-EDIT  SUBROUTINE  'GRAPH'). 

C 

C 

c . 

c 

C  ALPHABETICAL  LISTING  OF  VARIABLES 

C 

C  VARIABLE  (SUBROUTINE  USED)  DESCRIPTION. 

C 

C . 

c 

C  DELTAT . ( STEP , Q , OUTPUT , FORCE , GRAPH ) INCREMENTAL  STEP 

C  IN  NON-DIMENSIONAL  TIME  (INVERSE  OF  STEPS  PER 

C  CYCLE). 

C 

C  DFDZ . (QA,QB) COMPLEX  PARTIAL  DERIVATIVE  OF  THE  COMPLEX 

C  FUNCTION  F  WITH  TO  RESPECT  Z. 

C 

C  I . (QA,QB, FORCE) CONSTANT  COMPLEX  NUMBER  =  (0.0, 1.0) 

C 

C  ISTEPS . (ALL 'NUMBER  OF  INCREMENTAL  STEPS  PER  CYCLE  TO  BE 

C  TAKEN  FOR  ALL  ITERATIONS. 

C 

C  Kl-4 . (RUNGE)  RUNGE-KUTTA  CONSTANTS 

C 

C  KAPPAA . ( ALL ) STRENGTH  OF  THE  VORTEX  (GAMMA/ (2*PI*UMAX) ) 

n 

v- 

C  KAPPAB . (ALL) STRENGTH  OF  THE  VORTEX  (GAMMA/ ( 2* PI* UMAX) ) 
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.(ALL)  KEULLEGAN- CARP ENTER  NUMBER  FOR  THE  OSCILLATING 
FLOW. 


LI . (INPUT)  LOGICAL  OPERATOR. 


NUMCYC 


(INPUT) TOTAL  NUMBER  OF  CYCLES  TO  BE  EVALUATED. 


NUMPTS . (ALL)TOTAL  NUMBER  OF  POINTS  TO  BE  DETERMINED. 

(NUMCYC  *  ISTEPS) 


(ALL) CONSTANT  =  3.141592654 


QA . (QA) FUNCTION  SUBROUTINE  WHERE  Q  IS  THE  COMPLEX 

VELOCITY  OF  THE  VORTEX. 


QVECTA(N) . . 


.  (QB) FUNCTION  SUBROUTINE  WHERE  Q  IS  THE  COMPLEX 
VELOCITY  OF  THE  VORTEX. 

.  (STEP, GRAPH) COMPLEX  ARRAY  CONTAINING  THE  VELOCITY 
OF  THE  VORTEX  AT  ANY  POINT  Z(N). 


QVECTA(N) .... (STEP, GRAPH) COMPLEX  ARRAY  CONTAINING  THE  VELOCITY 
OF  THE  VORTEX  AT  ANY  POINT  Z(N) . 


TE . (RUNGE,QA,QB) NON -DIMENSIONAL  ENVIRONMENTAL  FLOW  TIME 


.  (STEP , QA , QB ) NON-DIMENSIONAL  TIME  CHANGED  TO 
RADIANS  FOR  COS/SIN  EVALUATIONS  (0-2~PI  FOR  A  CYCLE) 


TTF . (OUTPUT ) NON- DIMENSIONAL  FINISH  TIME  OF  THE  RUNS. 


.  (ALL) INITIAL  NON-DIMENSIONAL  STARTING  TIME  FOR  THE 
OSCILLATING  FLOW. 


.*  V  •  ‘  ‘  ‘ 


C  2A . (ALL)  INSTANTANEOUS  COMPLEX  LOCATION  OF  THE  VORTEX. 

C 

C  ZB . (ALL)  INSTANTANEOUS  COMPLEX  LOCATION  OF  THE  VORTEX. 

C 

C  ZAI . (GRAPH) IMAGINARY(Y)  LOCATIONOF  THE  VORTEX. 

C 

C  ZBI . (GRAPH)  IKAGINARY(Y)  LOCATIONOF  THE  VORTEX. 

C 

C  ZP . (QA,Q3) COMPLEX  PREDICTED  VALUE  FOR  Z(N) 

C 

C  ZAR . (GRAPH) REAL (X)  LOCATION  OF  THE  VORTEX. 

C 

C  ZBR . ( GRAPH) REAL (X)  LOCATION  OF  THE  VORTEX. 

C  ZZ . (QA , QB) DUMMY  ARGUMENT  FOR  THE  FUNCTION  SUBROUTINE  QA 

C  AMD  QB  WHICH  CONTAINS  THE  PASSED  VALUE  OF  Z(N). 

C  1 RUNGE )  DUMMY  ARGUMENT  CALL  FUNCTION  SUSROUITNE  QA 

CAND  QB 

C  . 

c 

c 

c 


C  *  MAIN  PROGRAM  * 

Q-kyr-rx+rk-kx-k-k  *rx  ■k'k-kjr'k-k-k-x-k-k-k  *  :r  rc  *  *  *tt  Ax*x^x  *******  **A****7x**:**3XX**:rxc**-**;A;7cx 

c 

IMPLICIT  REAL*8  (A-H,0-Z) 

REALMS  KC , KAFPAA , KAPPAB 
C0MPLEX*16  ZA ( 3649 ) , QVECTA ( 3649 ) 

COMPLEX* 16  Z3 ( 8649 ) , QVECTB ( 3649 ) 

COMMON  /FIRST/  TTI.ISTEPS  /SECOND/  NUMPTS  /THIRD/  KC, KAFPAA, 

1  KAPPAB /FOURTH/  PI 

c  \ 


C  SET  CONSTANTS 


KAPPAA  =  0.5D0 
KC  =  10. DO 
PI  =  3 . 141592654D0 
ISTEPS  =  1440 
C 

10  CALL  INPUT  (ZA , ZB ) 

CALL  STEP  (ZA,QVECTA,  ZB,QVECTB) 

CALL  OUTPUT ( ZA , ZB ) 

CALL  GRAPH  (ZA , QVECTA , Z3 , QVECT3 ) 

C 

WRITE  (6,*) 'DO  YOU  WISH  TO  RE-RUN  THE  PROGRAM?  1-YES,  O-NO1 
READ  (5,*)  LI 
IF  (Ll.NE.l)  GOTO  20 

WRITE  ( 5 , - )  'EXECUTION  CONTINUES  .  1 

GOTO  10 
20  CONTINUE 
C 

CALL  DONEPL 
STOP 
END 
C 

c 

C  *  SUBROUTINE  INPUT  * 

C  A**********************  -k-k  *******************  **************  A********** 

c 

C  THIS  SUBROUTINE  INTERACTIVELY  INPUTS  THE  VALUES  FOR  THE  INITIAL 
C  LOCATION  CF  THE  VORTEX,  STARTING  TINE,  ITERATION  STEP  SIZE,  AND 
C  THE  NUMBER  OF  CYCLES  TO  BE  EXAMINED. 

C 

SUBROUTINE  INPUT  (ZA,ZB) 

IMPLICIT  REAL*3  (A-H,0-Z) 

REALMS  KC,  KAPPAA,  KAPPAB 
COMPLEX* 16  ZA( 6649), ZB (3649) 

COMMON  /FIRST/  TTI , ISTEPS  /SECOND/  NUMPTS/THIRD/XC , KAPPAA , KAPPAB 


c 

C  READ  LAST  4  POINTS  OF  EACH  ARRAY  FOR  CONTINUATION  OPTIONS 
C 

REWIND  35 
DO  10  N  =  1,4 
10  READ (3 5,*)  ZA(N) ,ZB(N) 

READ (35 , *)  TTI 
READ (35 , *)  NUMCYC 
C 

C  CONTINUATION  OPTIONS 
C 

WRITE (6 ,  *)  'DO  YOU  WISH  TO  CONTINUE  THIS  RUN  FORM  THE  LAST? 
#  O-NO1 
READ (5,*)  LI 
IF  (L1.NE.0)  GOTO  20 
C 

WRITE (6, *) 1  INPUT  THE  VALUE  OF  X  :  ZA(1)=(X,X)' 

READ ( 5 , *)  X 

ZA ( 1 )  =  DCMPLX(X,X) 

C 

WRITE ( 6 ,  *)  'INPUT  THE  VALUE  OF  EPSILON  RADIUS:' 

READ (5,*)  ER 

X  =  ER  *  X 

ZB ( 1 )  =  DCMPLX (X , -X) 

C 

WRITE (6,*)  'INPUT  THE  VALUE  OF  EPSILON  KAPPA:' 

READ (5,*)  EK 
KAPPAB  =  KAPPAA  *  EK 
C 

TTI  =  0.5D0 
C 

C  CALL  SUBROUTINE  RUNGE  TO  GENERATE  FIRST  4  POINTS  FOR  EACH  ARRAY 
C 


1-YES 


CALL  RUNGE  (ZA,ZB) 


C  INPUT  THE  NUMBER  OF  CYCLES  TO  BE  ANALYZED 
C 


20  WRITE (6,*) 'THE  LAST  VALUE  FOR  NUMCYC  =',MUMCYC 

WRITE (6 ,*) 1  DO  YOU  WISH  TO  CHANGE  IT?  I-YES  O-NO' 

READ (5,*)  LI 
IF (LI. ME. I)  GOTO  30 

WRITE (6,*)  'INPUT  THE  NEW  VALUE  OF  NUMCYC : ' 

READ  (5,*)  N'JMCYC 
30  CONTINUE 
C 

NUMPTS  =  NUMCYC  *  I  STEPS  +  4 
C 

RETURN 

END 

C 

C 

Q  ******7r*yKjK******x****x**7r7rjt*jr*7'*************+.*K***  +  ir**-k-k7'-k-k7'+,+:-k-k-f'.'k  +  ±-k 

C  *  SUBROUTINE  RUNGE  * 

Q  iK******************************************  *******  ***7!  ****■**+:**■*: -XT*** 

c 

C  THIS  SUBROUTINE  USES  A  4TH  ORDER  RUNGE-KUTTA  METHOD  TO  GENERATE  THE 
C  FIRST  FOUR  VALUES  OF  Z. 

C 

SUBROUTINE  RUNGE (ZA,ZE) 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMPLEXn&  ZA  ( 8649 ),  Z3  (8649)  ,  QA ,  K1 ,  K2  ,  K3  ,  X4  ,  ZP1 ,  ZZ 
COMMON  /FIRST/  TTI,ISTEPS 
C 

DELTAT  =  1 .DO/DFLOAT(ISTEPS) 

C 


K1  =  DLETAT  *  QA(2Z , ZB (N) ,TE1 ) 

ZZ  =  ZA(N)+0 . 5D0*K1 

TE2  =  TE1+  0 . 5D0*DELTAT 

K2  =  DELTAT  *  QA(ZZ,ZB(N) ,TE2) 

ZZ  =  ZA(N)+0.5D0*K2 

K3  =  DELTAT  *  QA(ZZ,ZB(N),  TE2 ) 

ZZ  =  ZA(N)  +  K3 

TE4  =  TE1  +  DELTAT 

K4  =  DELTAT  *  QA(ZZ,ZB(N),  TE4) 

C 

ZA(N+1)  =  ZA(N)  +  (Kl-2 . D0*K2+2 . D0*K3+K4 ) /6 . DO 
C 

ZZ  =  ZB(M) 

K1  =  DLETAT  *  QB(ZA(N) ,ZZ,TE1) 

ZZ  =  ZB(N)+0 . 5D0*K1 

K2  =  DELTAT  *  QB (ZA(N) , ZZ , TE2 ) 

ZZ  =  ZB(N)+0.5D0*K2 

K3  =  DELTAT  *  QB(ZA(N),ZZ,  TE2) 

ZZ  =  ZB (N)  +  K3 

TE4  =  TE1  +  DELTAT 

K4  =  DELTAT  *  QB(ZA(M) ,ZZ,TE4) 

C 

5  ZB (N+l )  =  ZB(N)  +  (Kl-2 . D0*K2+2 . D0*K3+K4) /6 . DO 
C 
C 

RETURN 

END 

C 

C 

C  *  SUBROUTINE  STEP  * 

c 

C  THIS  SUBROUTINE  PROGRESSIVELY  CALCULATES  THE  VALUES  OF  ZA(N)  AND 
C  ZB (N)  EASED  ON  PREDICTOR-CORRECTOR  CONVECTION  SCHEME. 


.  • .  ■ .  »  •  ■'•y-y- »*»  .  •  ,n  *■» 
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SUBROUTINE  STEP  (ZA , QVECTA , ZB , QVECTB ) 

IMPLICIT  REAL*8  (A-H.O-Z) 

REALMS  KC , KAPPAA , KAPPAB 

COMPLEX* 16  ZA(8649 ) ,QVECTA(8649) ,ZP,QA,QB 
COMPLEX*16  ZB ( 8649 ), QVECTB (8649) 

COMMON  /FIRST/  TTI,ISTEPS  /SECOND/  NUMPTS/THIRD/  KC , KAPPAA , KAPPAB 
DELTAT  =  1  .DO/DFLOAT(ISTEPS) 

MAKE  VELOCITY  EVALUATIONS  FOR  POINTS  1  THROUGH  4 

DO  7  N  =  1,4 

TE  =  (DFLOAT  (N-l)  *  DELTAT  +  HI) 

QVECTA (N)  =  Q A ( ZA ( N ) , ZB ( N ) , TE ) 

7  QVECTB (N)  =  QB (ZA (N) , ZB (N) , TE ) 

DO  10  N  =  5,  NUMPTS 

NM1  =  N-l 

M2  =  N-2 

NM3  =  N-3 

NM4  =  N-4 

TE  =  ( DFLOAT (N-1)*DELTAT+TTI) 

CALCULATE  PREDICTED  VALUE  OF  ZA(N)  BASED  ON  ZB(N-l) 

Z?  =  ZA(NM1 )+( 55 . DQ*QVECTA(NM1 )  -59 . DO*QVECTA(NM2 ) 

1  +37.D0*QVECTA(NM3)-9.D0*QVECTA(HH4)) 

2  *DELTAT*  2. DO  *  KC  /  24. DO 

CALCULATE  1ST  ESTIMATE  OF  ZA(N)  BASED  ON  ZB(N-l) 

ZA(N)  =  ZA(MMl)  +  (9.D0*QA(ZP,ZB(NH1) ,TE)  +  19 . D0*QVICTA(NM1 ) 

1  -5.D0*QVECTA(NM2)+QVECTA(NM3)  ) *DELTAT  *  2. DO  *  KC  /  24. DO 


■-I''  -.  ZfZ?'.  ■■ Z ■  ■ 
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C  CALCULATE  PREDICTED  VALUE  OF  ZB(N)  BASED  ON  ESTIMATED  ZA(N) 
C 


ZP  =  ZB (NM1 ) + ( 55 . D0*QVECTB (NM1 )  -59 .D0*QVECTB (NM2 ) 

1  +37 . D0*QVECTB (NM3 ) -9 . D0*QVECTB (NM4) ) 

2  "DELTAT*  2. DO  *  KC  /  24. DO 
C 

C  CALCULATE  1ST  ESTIMATE  OF  ZB(N)  BASED  ON  ESTIMATED  ZA(N) 

Z3(N)  =  ZB(NMl)  +  (9.D0*QB(ZA(M),ZP,TE)  +  19.D0*QVECT3(NM1 ) 

1  - 5 . D0*QVECTB (NM2 ) +QVECTB (NM3 )  )*DELTAT  *  2. DO  *  KC  /  24. DO 

C 

C  CALCULATE  PREDICTED  VALUE  OF  ZA(N)  BASED  ON  ESTIMATED  ZB(N) 

C 

ZP  =  ZA ( NM1 }  +  ( 5 5 . DC *QVE CTA ( NM1 )  -59.D0*QVECTA(NM2) 

1  +37 . D0*QVECTA(NM3 ) -9 . D0*QVECTA(NM4) ) 

2  *DELTAT*  2. DO  *  KC  /  24. DO 
C 

C  CALCULATE  ACTUAL  ZA(N)  BASED  ON  ESTIMATED  ZB(N) 

C 

ZA(N)  =  ZA(NMl)  +  (9.D0*QA(ZP,ZB(N) ,TE)  +  19 . D0*QVECTA(NM1 ) 

1  -5.D0*QVECTA(NM2)+QVECTA(NM3)  )*DELTAT  *  2. DO  *  KC  /  24. DO 

C 

C  CALCULATE  PREDICTED  VALUE  OF  ZB(N)  BASED  ON  ACTUAL  ZA(N) 

C 

ZP  =  ZB(NMl)+(55. DO*QVECTB (NM1 )  -59 .D0*QVECTB(NM2) 

1  +37 . D0*QVECTB (NM3 ) -9 . D0*QVECTB (NM4 ) ) 

2  *DELTAT*  2. DO  *  KC  /  24. DO 
C 

C  CALCULATE  ACTUAL  ZB(N)  BASED  ON  ACTUAL  ZA(N) 

C 

ZB(M)  =  ZB(NMl)  +  ( 9 . D0*Q3 (ZA (N) , ZP , TE )  +  19.D0*QVECTB(NM1) 

1  - 5 . DQ*QVECTB ( NM2 ) +QVECTB (NM3 )  )*DELTAT  *  2. DO  *  KC  /  24. DO 

C 

C  EVALUATE  VELOCITIES  AT  ZA(N)  AND  ZB(N) 


QVECTA(N)  =  QA  ( ZA  ( N  )  ,  ZE  ( N  )  ,  TE  ) 
10  QVECTB(N)  =  QB (ZA (N) , ZB (N) ,TE) 


C 

RETURN 

END 

C 

C 


FUNCTION  SUBPROGRAM  QA 


C 

C  THIS  FUNCTION  SUBPORGRAM  IS  USED  TO  EVALUATE  VALUES  OF  QA 

C 

COMPLEX  FUNCTION  QA*16 ( ZA , ZB ,TE ) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL*3  KC,  KAPPAA , KAPPAB 
C0MPLEX*16  ZA , ZB ,  DFDZ ,  I,AA,AB,AC 

COMMON  /FIRST/  TTI,ISTEPS  /THIRD/  KC , KAPPAA , KAPPAB  /FOURTH/  PI 
C 

I  =  (0 .DO , 1 .DO) 

TT  =  TE  *  2. DO*  PI 
C 

AA  =  1 .DO/ (ZA-1 .DO/ (DCONJG(ZA) ) ) 

AB  =  1 .DO/ (ZA-ZB) 

AC  =  1 .DO/ (ZA-1 .DO/ (DCONJG(ZB) ) ) 

C 

DFDZ  =DSIN(TT)*(1 .DO-1 .DO/(ZA*ZA) ) -I*KAPPAA*AA-I*KA?PA3*( AB-AC) 

n 

QA  =  DCONJG(DFDZ) 

C 

RETURN 

END 

C 

r 
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c  I 

c  - 


FUNCTION  SUBPROGRAM  QB 


C 

C  THIS  FUNCTION  SUBPORGRAK  IS  USED  TO  EVALUATE  VALUES  OF  QB 
C 

COMPLEX  FUNCTION  Q3*16 (ZA, ZB , TE ) 

IMPLICIT  REAL* 3  (A-H,0-Z) 

REAL*3  KC,  KAPPAA , KAPPAB 
COMPLEX* 16  ZA , ZB ,  DFDZ ,  I,BA,BB,BC 

COMMON  /FIRST/  TTI,ISTEPS  /THIRD/  KC , KAPPAA , KAPPAB  /FOURTH/  PI 
C 

I  =  (0 . DO , 1 . DO ) 

TT  =  TE  *  2. DO*  PI 
C 

BA  =  1 .DO/ (ZB-ZA) 

BB  =  1 .DO/ (ZB-1 .DO/ (DCONJG(ZA) ) ) 

BC  =  l.DO/ (ZB-1. DO/ (DCONJG(ZB) ) ) 

C 

DFDZ  =DSIN(TT)*(1.D0-1.D0/(ZB*ZB))+I*KAPPAA*(BA-BB)+I*KAPPAB*BC 
C 

QB  =  DCONJG(DFDZ) 

C 

RETURN 
END 
C 
C 

C  *  SUBROUTINE  OUTPUT  * 

c 

C  THIS  SUBROUTINE  IS  USED  TO  SAVE  THE  FINAL.  VALUES  OF  Z.TT,  AND  NUNCYC 
C  FOR  POSSIBLE  CONTINUATION  OF  THE  PLOTS. 

C  THE  VALUES  ARE  SAVED  IN  FILE  'FILE  FT35F001 1 . 

C 

I 
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SUBROUTINE  OUTFUT(ZA,ZB) 

IMPLICIT  REAL*8  (A-H,  0-Z) 

COMPLEXES  ZA(8649) ,ZB(3649) 

COMMON  /FIRST/'  TTI,ISTEPS  /SECOND/  NUMPTS 
C 

DELTAT  =  1 .DO/FLOAT (ISTEP3) 

NN  =  NUMPTS  -  4 
TTF  =  D FLOAT (NN ) *  DELTAT  +  TTI 
NUMCYC  =  NN  /  I STEPS 
C 

REWIND  35 
DO  10  N  =  1,4 
NN  =  NN  +  1 

10  WRITE (35,*)  ZA(NN) , ZB(NN) 

WRITE (35,*)  TTF 
WRITE (35,*)  NUMCYC 
C 

RETURN 

END 

C 

c 

c 

c 

£  x**x*********x**xxx*x*******x***xx***xx******x*****x*******x*xx*x**** 

c  *  SUBROUTINE  GRAPH  * 

£  x*xx**xxxxx*xxxx*-**xx*xx***x*x*xxxxxxxx*x*xx*x*xxx*xxx**x*x****x*y l*i* 

c 

SUBROUTINE  GRAPH  -(ZA , QVECTA , Z3 , QVECTS ) 

IMPLICIT  REAL *3  (A-H,  O-Z) 

REAL*4  ZRA( 3649) ,ZIA( 8649) ,X( 361  ) , Y ( 3 6 1 ) , GARRAY ( 1 ) 

REAL *4  ZRB(6649) ,ZIB(8649) 

COMPLEX* 16  ZA ( 3649 ) , QVECTA (3649) 

COMPLEX* 16  ZB ( 8649 ) , QVECTB ( 8649 ) 

COMMON  /FIRST/  TTI,ISTEPS  /SECOND/  NUMPTS  /FOURTH/  PI 
C 
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C  DEFINE  CYLINDER  OF  RADIUS  1.0 
C 

DO  10  N  =  1,361 

RAD  =  FLOAT (N-l)  *  PI  / 180.0 

X(N)  =  COS (RAD ) 

10  Y (N)  =  SIM (RAD) 

C 

C  DEFINE  REAL  SINGLE  PRECISION  NUMBERS  FOR  PLOT  ROUTINES 
C 

DELTAT  =  1 . DO/DFLOAT ( I STEPS ) 

MPTSM4  =  NUMPTS  -  4 
C 

DO  20  N  =1,  NPTSM4 
ZRA(N)  =DREAL(ZA(N) ) 

ZIA(N)  =DIMAG(ZA(N) ) 

ZRB(M)  =DREAL(ZB(N) ) 

20  Z I 3 ( N )  =DIMAG(ZB(N) ) 

C 

C  GENERATE  PLOTS 
C 

C  CALL  COMPRS 

C  CALL  BLOWUP (1.0) 

CALL  TEK613 
CALL  ELGWUP (2.0) 

CALL  NOBRDR 
CALL  GRACE  (0.0) 

C 

C  VORTEX  TRACK  PLOT 

CALL  NOCHEK 

CALL  PAGE  (11.0,3.5) 

CALL  AREA2D ( 8 . 0 , 6 . 0 ) 

call  yna;:e(  '  ys  ,  100) 

CALL  XI LAKE  (  1  X3  '  ,  100) 

CALL  YAaAMG (0.0) 


CALL  INTAXS 


C 

CALL  XTICKS ( 2 ) 

CALL  YTICKS ( 2 ) 

CALL  GRAF (-26,1,2,-2,1,19) 

C 

CALL  CURVE ( X , Y , 3 6 1 , 0 ) 

C-ARRAY  ( 1 )  =  0.03 

CALL  SHADE (X,Y ,361,45, GARRAY ,1,0,0) 
C 

CALL  CURVE ( ZRA , 21 A , NPTSM4 , 0 ) 

C 

CALL  DASH 

CALL  CURVE ( ZR3 , Z IB , NPTSK4 , 0 ) 

CALL  EMDPL ( 0 ) 

C 

RETURN 

END 


APPENDIX  E 
FIGURES 


Figure  E.Ia  Experimental  In-line  Force  Trace 
Top  -  Force  Coefficient,  Bottom  -  Flow  Velocity 
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cure  E.5 


cure  E.6  Force  Coefficients  forZ(l)=(0, 
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Figure  E.l  1  u-Velocity  vs  x-Position  for  10  Cycles 
Ti  =  0.0 

(a)  Z(1)  =  (0,  1.50S3),  (b)  Z(I)  =  (0,  1.50S4) 


Y-POSITION 
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Figure  E.12  v-Velocity  vs  y-Position  for  10  Cycles 
Ti  =  0.0 


Figure  F..16  Velocity  Magnitudes  for  10  Cycles 
Ti  =  0  25 

(a)  7(1)=- (0.  1.6SS5).  (b)  Z(l)-(0,  1.6SS6) 
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Figure  E.  1 8 
(a)  Z(l) 


v-VcIocity  vs  y-Position  for  10  Cycles 
Ti  =  0.25 

=  (0,1.6885),  (b)  Z(1)  =  (0,  1.68S6) 


figure  E.27  Vortex  Path  Z(l)«(1.0.  1  0273) 
Ti  =  0  0 

(a)  Cvclcs  1-6,  (b)  Cycle  7,  (c)  Cvclcs  8-15 


Figure  E . 2 S  Velocity  Magnitudes  for  10  Cycles 
Ti  =  0.0 

(a)  Zm  =  U0.  1.02':),  (hi  Z(l)  =  (l  0.  1.0273) 
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